/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Definitions of the polychromatic_grid functions.

// $Id$

#include "polychromatic_grid.h"
#include "grid_factory.h"
#include "grid-fits.h"
#include "CCfits/CCfits"
#include "blitz-fits.h"
#include <vector>
#include "density_generator.h"
#include "arepo_grid.h"

using namespace CCfits;
using namespace std;
using namespace blitz;


template<template<typename> class T>
mcrx::polychromatic_grid<T>::polychromatic_grid (boost::shared_ptr<T_grid_impl> g) :
  g_(g)
{
  // can't call make_blocks here because there's no data in the
  // grid... it'll have to be called from derived constructors.

  //make_blocks ();
}

template<template<typename> class T>
mcrx::polychromatic_grid<T>::polychromatic_grid (boost::shared_ptr<T_grid_impl> g,
						 CCfits::ExtHDU& data_hdu,
						 const density_generator& gen,
						 int n_lambda) :
  g_(g)
{
  stupid(data_hdu, gen, n_lambda);
}


template<template<typename> class T>
void
mcrx::polychromatic_grid<T>::stupid (CCfits::ExtHDU& data_hdu,
				     const density_generator& gen,
				     int n_lambda)
{
  const int nc = n_cells ();
  cout << "Reading grid cell data for " << nc << " cells." << endl;

  Column& c_m_g = data_hdu.column("mass_gas");
  Column& c_p_g = data_hdu.column("p_gas");
  Column& c_m_metals = data_hdu.column("mass_metals");

  std::vector<T_float> m_g, m_metals;
  std::vector<vec3d> p_g;
  c_m_g.read(m_g, 1, c_m_g.rows() );
  c_m_metals.read(m_metals, 1, c_m_metals.rows() );
  read(c_p_g, p_g);

  // we need to convert momentum, which is in units of
  // mass*length/time, to units of mass*m/s (since what we want is the
  // velocity)
  const T_float momcon = units::convert(c_p_g.unit(),c_m_g.unit()+"*m/s");

  if (c_m_g.rows() != nc ||
      c_m_metals.rows() != nc) {
    cerr << "Error: HDU does not contain the expected number of rows" << endl;
  }
  assert(m_g.size()==nc);
  assert(m_metals.size()==nc);
  assert(p_g.size()==nc);

  int i = 0;
  for (iterator c = begin (), e=end(); c != e; ++c, ++i) {
    // make an absorber with the appropriate densities. We only know
    // density of metals and gas so we ask the density generator to
    // give us a density vector, which we put into our big block.
    assert(c->volume()>0);
    assert(m_g [i]>=0);
    assert(m_metals [i]>=0);
    T_densities rho (gen.make_density_vector(m_g [i]/c->volume(),
					     m_metals [i]/c->volume()));
    ASSERT_ALL(rho>=0);
    T_lambda intens (n_lambda);
    intens=0;
    T_data* d = new T_data (rho, p_g[i]*momcon/m_g[i], intens);
    DEBUG(1,cout << "Cell velocity " << d->velocity() << " in polychromatic_grid.cc" <<  endl;);
    c->set_data(d );
  }

  // and now consolidate into continuous blocks
  make_blocks();
}


/** This function goes through the grid and consolidates the cells and
    data to be in single blocks. (T_data copy constructor must be
    suitable for this.) */
template<template<typename> class T>
void mcrx::polychromatic_grid<T>::make_blocks()
{
  const int nc = n_cells ();
  // Allocate a block for our T_data, so they sit in a continuous
  // chunk of memory
  const size_t block_size = nc*sizeof (T_data);
  cout << "Allocating a memory block for " << nc 
       << " cell data objects, " << block_size*1.0/(1024*1024) << " MB."  << endl;
  data_block_ = operator new (block_size);

  // We also want to allocate blocks for the densities and intensities
  // arrays of the cells. For this we need to know the number of
  // species and wavelengths.
  const int n_species = begin()->data()->densities().size();
  const int n_lambda = begin()->data()->intensity().size();
  cout << "Allocating a (" << nc << ", " << n_species
       << ") memory block for cell densities, "
       << nc*n_species*sizeof(T_float)*1.0/(1024*1024*1024) << " GB.\n";
  rho_.resize(nc, n_species);
  cout << "Allocating a (" << nc << ", " << n_lambda
       << ") memory block for cell intensities, "
       << nc*n_lambda*sizeof(T_float)*1.0/(1024*1024*1024) << " GB.\n";
  intensities_.resize(nc, n_lambda);
  intensities_=0.;

  // go through and change the data objects
  int i = 0;
  T_data* block_pointer = reinterpret_cast<T_data*> (data_block_);
  for (iterator c = begin (), e=end(); c != e; ++c, ++i) {
    // and now we initialize the absorber with WEAK REFERENCES to the blocks
    // so that there is no reference counting them
    T_densities rho;
    rho.weakReference(rho_(i,Range::all()));
    T_lambda intens;
    if (n_lambda>0) 
      // if n_lambda==0, there is nothing to reference...
      intens.weakReference(intensities_(i,Range::all()));
    // and copy data from existing absorber
    rho = c->data()->densities();
    if(n_lambda>0) {
      intens.weakReference(intensities_(i,Range::all()));
      intens = c->data()->intensity();
    }
    ASSERT_ALL(intens==intens);
    ASSERT_ALL(intens< 1e300);
    // use placement new to make a copy of the data object 
    T_data* d = new (block_pointer++) T_data (rho, c->data()->velocity(), 
					      intens);
    c->set_data(d );
  }
  ASSERT_ALL(intensities_==intensities_);
  ASSERT_ALL(intensities_< 1e300);

}



/** Destructor.  The destructor must check if we have used placement
    new to allocate the cells and explicitly call the destructors for
    all the cell data if that's the case.  (Note that the destructor
    for the mcrx::adaptive_grid class does the same for the subgrids
    and the grid cells, but the cell data is stored separately.) */
template<template<typename> class T>
mcrx::polychromatic_grid<T>::~polychromatic_grid ()
{
  if (data_block_) {
    // The data was allocated with placement new, so we can't let the
    // grid delete them when it exits.  Explicitly call the
    // destructors for them now.
    for (iterator c = begin (); c != end(); ++c) {
      c->data()->~T_data();
      c->set_data_nodelete (0 );
    }
    operator delete (data_block_);
  }

  // The array blocks will selfdestruct implicitly and this is safe
  // because at this point, the weak references in the absorber
  // objects have already been destroyed.
}


// Explicitly intantiate the templates necessary to avoid having to
// include everything
template class mcrx::polychromatic_grid<mcrx::adaptive_grid>;

#ifdef WITH_AREPO
template class mcrx::polychromatic_grid<mcrx::arepo_grid>;
#endif
